# Estimate the conditional effects of populist involvement using alternative cutoffs to define populist
# parties (populist score greater than 8, 8.5, or 9). See Figure N.1, Appendix N.




#############
# load data #
#############

library(brms)
library(ggplot2)
library(ggpubr)
library(stargazer)

load("PSRM Replication Files/PresidentialElectionsPosts.RData")

colnames(president_df)

all(c("populist_present80", "polarization", "populist_present85", "populist_present90") %in% colnames(president_df)) # TRUE

length(unique(president_df$country)) # 26
length(unique(president_df$electionname)) # 29
length(unique(president_df$party_country)) # 52

# focus on 15 days before and after the election
day_range <- 15

sub <- president_df[which(president_df$daysinceelection >= -1*day_range
                          & president_df$daysinceelection <= day_range),]

# winners and losers
winners <- sub[which(sub$win == 1),]
losers <- sub[which(sub$win != 1),]

length(unique(winners$electionname[which(!is.na(winners$polarization))])) # 27
length(unique(winners$electionname[which(!is.na(winners$economic_polarization))])) # 27
length(unique(winners$electionname[which(!is.na(winners$social_polarization))])) # 27
length(unique(winners$electionname[which(!is.na(winners$populist_present))])) # 27
length(unique(winners$electionname[which(!is.na(winners$populist_dummy))])) # 27

length(unique(losers$electionname[which(!is.na(losers$polarization))])) # 27
length(unique(losers$electionname[which(!is.na(losers$economic_polarization))])) # 27
length(unique(losers$electionname[which(!is.na(losers$social_polarization))])) # 27
length(unique(losers$electionname[which(!is.na(losers$populist_present))])) # 27
length(unique(losers$electionname[which(!is.na(losers$populist_dummy))])) # 27




##################
# winners + love #
##################

winner_love80 <- brm(loveprop ~ post*populist_present80 + post*polarization
                     + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                     + Xl + Xr 
                     + (1 + Xl + Xr | party_election),
                     data=winners,
                     warmup=1000, iter=2000, chains=3, seed=1234)
max(summary(winner_love80)[[14]][,5])
#plot(winner_love80)
#summary(winner_love80)

winner_love85 <- brm(loveprop ~ post*populist_present85 + post*polarization
                     + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                     + Xl + Xr 
                     + (1 + Xl + Xr | party_election),
                     data=winners,
                     warmup=1000, iter=2000, chains=3, seed=1234)
max(summary(winner_love85)[[14]][,5])
#plot(winner_love85)
#summary(winner_love85)

winner_love90 <- brm(loveprop ~ post*populist_present90 + post*polarization
                     + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                     + Xl + Xr 
                     + (1 + Xl + Xr | party_election),
                     data=winners,
                     warmup=1000, iter=2000, chains=3, seed=1234)
max(summary(winner_love90)[[14]][,5])
#plot(winner_love90)
#summary(winner_love90)




###################
# winners + angry #
###################

winner_angry80 <- brm(angryprop ~ post*populist_present80 + post*polarization
                      + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                      + Xl + Xr 
                      + (1 + Xl + Xr | party_election),
                      data=winners,
                      warmup=1000, iter=2000, chains=3, seed=1234)
max(summary(winner_angry80)[[14]][,5])
#plot(winner_angry80)
#summary(winner_angry80)

winner_angry85 <- brm(angryprop ~ post*populist_present85 + post*polarization
                      + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                      + Xl + Xr 
                      + (1 + Xl + Xr | party_election),
                      data=winners,
                      warmup=1000, iter=2000, chains=3, seed=1234)
max(summary(winner_angry85)[[14]][,5])
#plot(winner_angry85)
#summary(winner_angry85)

winner_angry90 <- brm(angryprop ~ post*populist_present90 + post*polarization
                      + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                      + Xl + Xr 
                      + (1 + Xl + Xr | party_election),
                      data=winners,
                      warmup=1000, iter=2000, chains=3, seed=1234)
max(summary(winner_angry90)[[14]][,5])
#plot(winner_angry90)
#summary(winner_angry90)




#################
# losers + love #
#################

loser_love80 <- brm(loveprop ~ post*populist_present80 + post*polarization
                   + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                   + Xl + Xr 
                   + (1 + Xl + Xr | party_election),
                   data=losers,
                   warmup=1000, iter=2000, chains=3, seed=1234)
max(summary(loser_love80)[[14]][,5])
#plot(loser_love80)
#summary(loser_love80)

loser_love85 <- brm(loveprop ~ post*populist_present85 + post*polarization
                    + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                    + Xl + Xr 
                    + (1 + Xl + Xr | party_election),
                    data=losers,
                    warmup=1000, iter=2000, chains=3, seed=1234)
max(summary(loser_love85)[[14]][,5])
#plot(loser_love85)
#summary(loser_love85)

loser_love90 <- brm(loveprop ~ post*populist_present90 + post*polarization
                    + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                    + Xl + Xr 
                    + (1 + Xl + Xr | party_election),
                    data=losers,
                    warmup=1000, iter=2000, chains=3, seed=1234)
max(summary(loser_love90)[[14]][,5])
#plot(loser_love90)
#summary(loser_love90)




##################
# losers + angry #
##################

loser_angry80 <- brm(angryprop ~ post*populist_present80 + post*polarization
                     + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                     + Xl + Xr 
                     + (1 + Xl + Xr | party_election),
                     data=losers,
                     warmup=1000, iter=2000, chains=3, seed=1234)
max(summary(loser_angry80)[[14]][,5])
#plot(loser_angry80)
#summary(loser_angry80)

loser_angry85 <- brm(angryprop ~ post*populist_present85 + post*polarization
                     + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                     + Xl + Xr 
                     + (1 + Xl + Xr | party_election),
                     data=losers,
                     warmup=1000, iter=2000, chains=3, seed=1234)
max(summary(loser_angry85)[[14]][,5])
#plot(loser_angry85)
#summary(loser_angry85)

loser_angry90 <- brm(angryprop ~ post*populist_present90 + post*polarization
                     + incumbentparty + concurrent + runoff + semipresidential + enc_v1
                     + Xl + Xr 
                     + (1 + Xl + Xr | party_election),
                     data=losers,
                     warmup=1000, iter=2000, chains=3, seed=1234)
max(summary(loser_angry90)[[14]][,5])
#plot(loser_angry90)
#summary(loser_angry90)

#save(winner_love80, winner_love85, winner_love90,
#     winner_angry80, winner_angry85, winner_angry90,
#     loser_love80, loser_love85, loser_love90,
#     loser_angry80, loser_angry85, loser_angry90,
#     file="Stan_PopulismAlternativeCutoffs.RData")




##################
# result summary #
##################

load("PSRM Replication Files/Stan_PopulismAlternativeCutoffs.RData")

# plot interaction terms
getStdPostCI <- function(model, var1, var2, intr, label){
  coeftab <- summary(model)[[14]][,c(1,3,4)]
  coeftab <- coeftab[intr,]
  
  temp <- data.frame(label=label,
                     est=coeftab[1],
                     lwr=coeftab[2],
                     upr=coeftab[3])
  colnames(temp) <- c("label", "est", "lwr", "upr")
  return(temp)
}

citab <- rbind(getStdPostCI(winner_love80, 2, 3, 12, "Winner + Love"),
               getStdPostCI(winner_love80, 2, 4, 13, "Winner + Love"),
               getStdPostCI(winner_love85, 2, 3, 12, "Winner + Love"),
               getStdPostCI(winner_love85, 2, 4, 13, "Winner + Love"),
               getStdPostCI(winner_love90, 2, 3, 12, "Winner + Love"),
               getStdPostCI(winner_love90, 2, 4, 13, "Winner + Love"),
               getStdPostCI(winner_angry80, 2, 3, 12, "Winner + Angry"),
               getStdPostCI(winner_angry80, 2, 4, 13, "Winner + Angry"),
               getStdPostCI(winner_angry85, 2, 3, 12, "Winner + Angry"),
               getStdPostCI(winner_angry85, 2, 4, 13, "Winner + Angry"),
               getStdPostCI(winner_angry90, 2, 3, 12, "Winner + Angry"),
               getStdPostCI(winner_angry90, 2, 4, 13, "Winner + Angry"),
               getStdPostCI(loser_love80, 2, 3, 12, "Loser + Love"),
               getStdPostCI(loser_love80, 2, 4, 13, "Loser + Love"),
               getStdPostCI(loser_love85, 2, 3, 12, "Loser + Love"),
               getStdPostCI(loser_love85, 2, 4, 13, "Loser + Love"),
               getStdPostCI(loser_love90, 2, 3, 12, "Loser + Love"),
               getStdPostCI(loser_love90, 2, 4, 13, "Loser + Love"),
               getStdPostCI(loser_angry80, 2, 3, 12, "Loser + Angry"),
               getStdPostCI(loser_angry80, 2, 4, 13, "Loser + Angry"),
               getStdPostCI(loser_angry85, 2, 3, 12, "Loser + Angry"),
               getStdPostCI(loser_angry85, 2, 4, 13, "Loser + Angry"),
               getStdPostCI(loser_angry90, 2, 3, 12, "Loser + Angry"),
               getStdPostCI(loser_angry90, 2, 4, 13, "Loser + Angry"))

citab$cutoff <- as.factor(rep(c(8, 8.5, 9), each = 2))

citab$intr <- c("(a) Post Election × Populist Involvement", "(b) Post Election × Polarization")

citab$label <- factor(citab$label, levels=rev(unique(citab$label)))
citab$intr <- factor(citab$intr, levels=unique(citab$intr))

# figure N.1
ggplot(data = citab, 
       aes(x = est, y = label, group = cutoff, color = cutoff)) + 
  facet_wrap(. ~ intr, ncol = 2, scales = "free_x") +
  geom_errorbar(data = citab, aes(xmin = lwr, xmax = upr),
                lwd=0.8, width = 0.2, position = position_dodge(width = -0.4)) +
  geom_vline(xintercept=0, lwd=0.5, lty=2) +
  geom_point(cex=4, pch=19, position = position_dodge(width = -0.4)) +
  scale_color_grey(name = "Populism Cutoff") +
  xlab("Posterior Mean and 95% Credible Interval of Interaction Term") +
  ylab("") +
  theme_bw() +
  theme(strip.text.x = element_text(size = 14),
        axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12),
        axis.title.x = element_text(size = 16),
        axis.title.y = element_text(size = 16),
        panel.spacing = unit(1, "lines"))

#ggsave("interaction_populismcutoff.pdf", width=12, height=7)
